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Abstract 

Wc consider a self-propcllcd particle system which has been used to describe cer- 
tain types of collective motion of animals, such as fish schools and bird flocks. 
Interactions between particles are specified by means of a pairwise potential, 
repulsive at short ranges and attractive at longer ranges. The exponentially 
decaying Morse potential is a typical choice, and is known to reproduce cer- 
tain types of collective motion observed in nature, particularly aligned flocks 
and rotating mills. We introduce a class of interaction potentials, that we call 
Quasi-Morse, for which flock and rotating mills states are also observed nu- 
merically, however in that case the corresponding macroscopic equations allow 
for explicit solutions in terms of special functions, with coefficients that can be 
obtained numerically without solving the particle evolution. We compare the 
obtained solutions with long-time dynamics of the particle systems and find a 
close agreement for several types of flock and mill solutions. 
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1. Introduction 

Emerging behaviors in interacting particle systems have received a lot of 
attention in research in recent years. Topics range from diverse fields of appli- 
cations such as animal collective behavior, traffic, crowd dynamics and crystal- 
lization. Self-organization in the absence of leaders has been reported in several 
species which coordinate their movement (swarming) , and several models have 
been proposed for their explanation 



Many of these models are based on zones in which some of 3 basic effects are 
included: short-range repulsion, long-range attraction, and alignment. These 
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3-zone basic descriptions have been very popular for modeling fish schools [30|, 
32, H,!^, starlings [2^, or ducks [li, 35|. The main modelling issues are if 



some or all of these effects between agents have to be included and if so, how 
to incorporate them. Many basic swarming models rely on averaged spatial 
distance or orientation interactions while recent biological studies point out the 
importance of nearest-neighbor interactions 0] or anisotropic communication 
[31 1 . Mathematicians have started in recent years to attack one of the most 
striking features of these simple looking models: the diversity of swarming states, 
also called patterns in the biology community, their emergence and stability. 

The individual level description of these phenomena leads to certain particle 
systems, called Individual Based Models (IB Ms), with some common aspects. 
Typically, the attraction-repulsion is modeled through pairwise efi^ective poten- 
tials depending on the distance between individuals. An asymptotic speed for 
particles is imposed either by working in the constrained set of a sphere in ve- 
locity space [U [25I, [1^ or by adding a term of balance between self-propulsion 
and friction which effectively fixes the speed to a limiting value for large times 
In this work, we will not include any alignment mechanism. Wc refer 



to jl3| for a survey on results related to kinetic modeling in swarming. 

In Section 2 we will review some of these IBMs, and discuss the appear- 
ance of two main swarming patterns: mills and flocks. These patterns are easily 
observed in particle simulations [II, 11| and reported in detail for certain partic- 
ular potentials, the so-called Morse potentials. We will give a precise definition 
of flocks and mills as solutions of the kinetic equation associated to the particle 
systems. Finding the spatial shape of flocks and mills has been numerically re- 
ported in the literature but obtaining analytical results on them has only been 
done in one dimension for the Morse potential in . 

In this work, we generalize the strategy in [HI proposing a new interaction 
potential, that wc call Quasi-Morsc, to replace the Morse potential. The Quasi- 
Morse potential coincides with the Morse potential in one dimension and wc will 
show that it is a suitable extension of the Morse potential in 71 = 2, 3. Section 
3 introduces Quasi-Morse potentials as fundamental solutions of certain linear 
PDEs. We will first show that the Quasi-Morse potentials are biologically rele- 
vant in essentially the same parameter range as the Morse potentials. Second, 
we make use of their particular structure to show in our main theorem that 
flock and mill solutions can be expressed as almost explicit linear combinations 
of special functions. 

Finally, Section 4 is devoted to propose an algorithm to compute the scalar 
coefficients in the expansion of the flock and mill patterns in terms of the basis 
functions associated with the Quasi-Morse PDE operators. The strategy uses 
ideas of constrained optimization methods. We finally compare the results for 
flocks in 2D and 3D and mills in 2D to particle simulations showing a good 
agreement. As a conclusion, we demonstrate that the proposed Quasi-Morse 
potentials are a very good alternative to Morse potentials as they share many 
of their features in the natural parameter range, and at the same time enable 
explicit computation of the macroscopic density profiles up to numerically de- 
termined constants. 
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2. Swarming: Models & Patterns 



We will consider a simple second order model for swarming analyzed in [21 1 
consisting of the attraction-repulsion of N interacting self-propelled particles 
located at Xi S R" with velocities Vi S R" in a host medium with friction, 
with n = 1,2,3. Friction is modeled by Rayleigh's law and as a result, an 
asymptotic speed for the individuals is fixed by the compensation of friction and 
self-propulsion. More precisely, the time evolution is governed by the equations 
of motion 



~dt 

dv. 

~dt 



= avi - I3vi\vi\^ ~ ^ W{xi - Xj) , 



(1) 



potential and a, /3 are effective values for 
16, [l^ for more discussion. The 
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where is a pairwise interaction 
propulsion and friction forces, see [33 
interaction potential W : R" x — > M is assumed to be radially symmetric: 
W{x) = U{\x\), X e M". The typical asymptotic speed of the individuals is 
The Morse potential is defined by taking 



U{r) = -C^e-''/'-* 4- C_ 



Re 



-r/lR 



where Ca, Cr are the attractive and repulsive strengths, and Ia, Ir are their 
respective length scales. We set V{r) = — exp(— r/Z^i), C = Cr/Ca, and / = 
Ir/Ia to obtain 



U{r) = Ca [Vir) - CV (y 



The choice of this potential is motivated in [21| for being one of the simplest 
choices of integrable potentials with easily computable conditions to distinguish 
the relevant parameters in biological swarms. In fact, it is straightforward to 
check that in the range C > 1 and / < 1 the potential U{r) is short-range 
repulsive and long-range attractive with a unique minimum defining a typical 
distance between particles. Moreover, in this regime the sign of the integral of 
the potential: 



w{x) dx ^ v(i - cr) 



with V := 



V^(r)r"-idr < 0, (2) 



gives a criterion to distinguish between the so-called H-stablc and catastrophic 
regimes. This condition reads as CI" — 1 < for the catastrophic case in any 
dimension n. 



see 
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39| . This property of the potential is important since it is 
related to the typical patterns emerging in such systems, as classified in (2]| . 

Flocks, where particles tend to form groups, moving with the same velocity, 
and milling solutions, where rotatory states are formed are of par ticular interest 
and are observed in particle and hydrodynamic simulations |2ll . 



14 in n = 2. 



Actually, they typically emerge in the large time behavior of the system of 
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particles ([T]) in the catastrophic regime CP < 1 with C > 1 and ^ < 1. In 
the same range of parameters, randomly chosen initial data lead also to other 



patterns such as double mills and flocks [2l|, [12| . However mills are not observed 
in the H-stable regime CP > 1 with C > 1 and I < 1 while flocks do. 

Assuming the weak coupling scaling [l^, [s^, H, in which the range of 
interaction is kept fixed and the strength of interaction is divided proportionally 
between particles, we pass to the rescalcd formulation: 

dxi 



This system has a well-defined limit as iV — > oo which can be expressed as a 
solution of the corresponding mean-field equation: 

dtf + v- V,/ + F[p] ■ V„/ + div ((a - P\v\^) «/) = , (3) 

with 

p{t,x) := / f{t,x,v)dv. 



Here, f{t,x,v) : R x M" x M" ^ M is the phase-space density, and p{t,x) 
is the averaged (macroscopic) density. The mean-field interaction is given by 

The limit N ^ oo has been established rigorously for smooth potentials 
W e in [iO, m, H 113], in [l^, Q for more general models with and without 



noise, and for more general potentials, with possibly singular behavior at zero. 



including the Morse potential ([2]) in the recent result 27|. 



2.1. Flock and Mill States 

We are interested in computing certain relevant particular solutions of the 
Vlasov-like equation for swarming in (|3]). In fact, we can formally find mono- 
kinetic solutions of ([3]) by inserting the ansatz: 

f{t,x,v) = p{t,x)6{v -u{t,x)), (4) 

in the weak formulation of (jSj. The result in [3, is that p and u should 
satisfy the following set of hydrodynamic equations: 

dp 

— + div^{pu) = 0, 

(5) 

P{u- "^x)u = p{a- (3\u\^)u - p (V^W ★ p). 

ot 

Definition 1. A flock is a solution fp of ([3]) of the form ^ with p(t,x) — 
Pf{x — tuo) and u{t,x) — uq with uq € M" such that \uo\ — and pp a 

probability measure in R". 
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Obviously, flock solutions are determined by their density profile pp and 
have the structure of traveling waves in the direction of the velocity vector uq. 
It is straightforward to sec that the density of a flock is characterized by the 
following equation: 

Proposition 1. The Junction fp > is a flock solution if and only if the 
macroscopic density pp satisfies 

-k pp ~ on the support of pp . (6) 

There arc singular solutions to (jHl) obtained by concentrating all the mass uni- 



formly in a suitable sphere, the so-called Delta rings [12|, whose stability for 
first order models has recently been studied in [l[ for certain potentials. Also, 
there are solutions to ([6|) given by smooth compactly supported densities for 
combination of suitable powers in ID [22|, [2^, for the Morse potential in ID 
, and for combination of powers when one of them is the repulsive Newtonian 



potential |2J] in 2D. In fact, the set of solutions to ^ can be very complicated 



even in one dimension [2^ |23|, [26| depending on the regularity of the potential. 

Let us remark that since we assume the radial symmetry of the potential, one 
expects that the density of the flocking solutions to ^ is radially symmetric as 
well and that it is supported in a ball B{0, Rp) with Rp > 0. This is reinforced 
by the fact that the convolution of radial functions is radial, see Subsection 2.2 
for more precise statements. Wc will reduce ourselves to find flocking solutions 
with radial symmetry in the rest of this paper, that is, finding Rp > and a 
radial density compactly supported in B(0,Rp) satisfying 

W*pp^C inB(0,i?F), (7) 

for some constant C € M. 

Another interesting type of solutions that spontaneously show up in particle 
simulations arc mills, they correspond to motion with the velocity field of a 
point vortex: 

.M(.) = ±yf:^, (8) 

where x = {xi,X2), x-^ ~ {~X2,xi), such that /9m(|2^|) is a radially symmetric 
stationary solution to ([5]). 

Definition 2. A mill is a solution fi\i of Q of the form: 

fM{t,x,v) = pm{x)5{v - um{x)) , 
with um given by (|8]) and pM radially symmetric. 

As shown in 0, solutions can also be characterized as: 

Proposition 2. Pm{x) is a mill density if and only if 



OL 

W-kp- -logl^l 



0, on the support of p . 
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As discussed above, one can obtain singular mill solutions by concentrating 
all particles in a ring [12[. However, we will search for radial solutions supported 
in an annulus B{Rm,RM) with < Rm < Rm, and therefore, mill radial 
solutions supported in B{Rrm Rm) are characterized by 

W^^PA/ = i^ + -log|a;| in B{Rm,RM), (9) 

for some constant D G R. In the following, the subindex x in differential 
operators is dropped since we only deal with x-dependent functions. 

2.2. Convolution of radial functions 

Since we want to find particular radial solutions to flocks ([7]) and mills ©, we 
need suitable expressions of the convolution of two radial functions in n = 2, 3. 
Given any radial density /3(|a;|), then the convolution term rewrites: 

{W-kp){x)= [ W{x~y)p{\y\)dy= [ [ W{x ~ SLj)p{s)s''~^dujds 

Jr" Jo JdB{0.1) 

which is not a convolution in r = |a:;| anymore, but rather is given by an integral 
operator of the following form: 



{W*p){r) = / 1'(r,s)p(s)ds 

Jr+ 

with 

*(r, s) = s"-i [ U{\rei - sa;|)da; . 

JdB{0.1) 

Expressing it in polar {n — 2) or spherical (n = 3) coordinates, we get the 
functions 

/•27r 

^{r,s) = sj U i^Vr^ -"^rs cos 9 + s^j dd (10) 

for n ~ 2 and 

*(r,s) = s^/ / U{\rei - suj{d,i^)\) smiydiyd9 
Jo Jo 

= 2tts^ J U (^\^r'^ — 2rs cos ly + s^^ siui/dz^, (11) 
with u!{0, v) — (cos V, sin v cos 9, sin v sin 0) for n = 3. 



3. Quasi-Morse potentials and their explicit solvability 

In this section, we define Quasi-Morse potentials for n = 1, 2, 3 and discuss 
their properties. These Quasi-Morse potentials will yield biologically relevant 
shapes similar to the Morse potentials. We show that flock and mill solutions 
in the natural parameter range, see Figures |4] and [IJd) for precise statements, 
can be computed explicitly up to numerically determined constants. 
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(a) Quasi-Morse potential 



(b) Morse potential 



Figure 1: Comparison of potentials: Both yield the biologically relevant shape 
of short-range repulsion and long-range attraction (Quasi-Morse: 71 = 2, C = 
fj = 0.75, fc = i, A = 4, Morse: C = f,l = 0.75, A = 2). 

3.1. Definition and comparison 

Definition 3. Let V : — > M denote the radially symmetric solution of the 
n-dimensional screened Poisson equation Au — k^u = Sq, for a given k > 0, 
that vanishes at infinity. Let C,l,X G M be further positive parameters. Then 
we say that U{\x\) is the n-dimensional Quasi-Morse potential if 

U{r) :^\(y{r)-Cv{^-)) . 

Using the radially symmetric ansatz, the screened Poisson equation reduces 
to a second-order ordinary differential equation dependent on the space dimen- 
sion. For relevant n = 1,2,3 this ODE possesses two linearly independent 
solutions. We therefore have 

Corollary 1. Quasi-Morse potentials for n = 1,2,3 are well-defined and con- 
structed from the following fundamental solution: 

71 = 1 : V{r) = -^e-^"" 

n^2: V{r) = -^Ko{kr) (12) 
n = 3: y(r) = -^^ 

where Kq is the modified Bessel function of second kind. For n = I, the Quasi- 
Morse potential equals the Morse potential. 

We illustrate the Quasi-Morse potential in comparison to the Morse potential 
for n = 2 with parameters C ~ 10/9 and / = 0.75 in Figure [1] Both potentials 
could be used to model the biologically motivated interplay between short-range 
repulsion and long-range attraction, and there is no clear reason to prefer one 
over the other. A significant difference is the behavior at zero, where Morse is 
finite and Quasi-Morse is singular though locally integrable for tt, > 1, which 
are the dimensions we aim to study. The parameter dependence of catastrophic 
regimes is inherited from the Morse potential as summarized in the next result, 
whose proof is given in an appendix. 
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Corollary 2. The function U{r) has a unique minimum if and only if I < 1, 
(jin-2 ^ ^ Furthermore, the Quasi-Morse potential U{\x\) is catastrophic if 

cr' < 1. 

Remark 1. We first emphasize that the global minimum of U corresponds to 
the biologically relevant scenario of short-range repulsion and long-range attrac- 
tion, as for the standard Morse potential. Concerning the H-stability of the 
Quasi-Morse potentials, we remark that the inverse Fourier transform of U (r) 
for k ~ 1 reads 

- ^ ci--i+p{ci"~^-m\^ 
^1^1^ (i + i^|2)(i+p|^i2) ■ 

which is positive if Cl"" > 1 and Cl""^^ > 1. This indicates the H-stability, 
but the criteria developed in fs^ l do not apply directly, since U{\S,\) is not inte- 
grable in dimensions rt = 2,3. However, our numerical findings presented in the 
following sections will suggest H-stability for the configurations I < 1,C/"^^ > 
1,C/" > 1. This corresponds to potentials having a unique minimum and a 
positive n— dimensional integral. 

Next, we mention the influence of the free scaling parameter k and show that 
any potential shape can be normalized to k = 1, see appendix. The following 
results are given without proof which follows easily by a change of variables 
from the convolution form in radial coordinates (jlOp and (fTT|) in subsection 2.2. 

Corollary 3. Let p be the flock (resp. mill) solution setting k ~ I with support 
-6(0, i?) (resp. B{Rm, Rm)), then the transformed solution p for the potential 
scaled to k = k ^ 1 is given by 

flock: p{x) = k'^p{kx) , supp(/5) = B(0, |) 
mill: p{x) ^ Pp{kx) , supp(/5) = B{S^, ^) ' 

Denoting U{r) := U{kr), W{x) = U{\x\), we have W -k p -k p = Ck"^'" for 
flocks, and W-kp = W-kp-\- log(fc) for mill solutions. 

3.2. Explicit solvability 

In this section, we show how to solve almost explicitly the integral equations 
for flock and mill profiles with the Quasi-Morse potential. The exact problem 
to solve for any potential W is to find a density p and its support such that 

(W * p){r) = s{r) on supp(p) (13) 

with some radial s(r) on supp(p) = B{Rjn, Rm), < Rm < Rm- Solving 
generally implies inverting the integral operator, a task that is complicated by 
the fact that the support is unknown. Even if this is numerically achievable, we 
will not learn anything about the structure of the solutions. For the Quasi-Morse 
potential, we take advantage of the differential operators behind its construc- 
tion, to avoid the inversion of (|13|) and to give an almost explicit expression of 
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its solution in terms of special functions. This strategy was already done in [5| 
in the one dimensional case, where Morse and Quasi-Morse potentials coincide. 
In this section, we pursue a similar strategy for the Quasi-Morse potential in 
n = 2,3. 

We begin our discussion reminding the radially symmetric fundamental sys- 
tem associated with some operators, which will be needed henceforth. 



Remark 2. The n- dimensional Helmholtz equation reads Au + /s^u = m R". 
Its fundamental system of radially symmetric solutions {ipi,ip2\ is associated to 
a second-order ordinary differential equation (in radial coordinates for n = 2,5) 
and given below, together with the fundamental system of the already mentioned 
screened Poisson equations {tpi,ip2}: 



Helmh. 




^1 


s.Poiss. 




1p2 


n = 1 




— ^ cos(fcr) 


n = l 




i p-kr 
2k 


n = 2 




^Yoikr) 


71 2 






n = 3 


1 sin(A;r) 


1 cos(fcr) 


n = 3 




1 e->"- 


4^ r 


47r r 


47r r 


47T r 



Here, Jq andYo are the Bessel functions of the first and second kind respectively, 
while Iq and Kq are the modified Bessel functions of the first and second kind 
respectively. 



We continue with a simple computation related to the local properties of our 
potential. 



Lemma 1. Let V be the fundamental solution of the screened Poisson equation. 
Then 

(''(7))4" + '""*- 



A 

Proof. Let ^ be a test function. Then by change of variables 



k^V{z)^{lz)dz + i{0) 



1 



R" 

fcV(z)^(?z)/"dz + r-2^(o) 



leading to the weak formulation of the claim. 
Now, we can state the main result of this section. 



9 



Theorem 1. Assume there exists a solution of {W * p){r) = s{r) on supp(p) 
with W being the Quasi-Morse potential and swpp(p) ~ B{0, Rp) , s(r) = D 
for flocks, or supp(yo) = B{Rm,RM) , s{r) = D + ^ log(r) for mills respectively. 
Then p has to be of the following form on supp p : 



n = 2; 


flock 


^ > 


Pf = Pi Jo{ar) + p2 






v4 = 


Pf = /ii?-^ + P2 






^ < 


Pf = Pi Io{ar) + P2 




mill 


^ > 


Pm = Pmhoni + Pi Joiar) + p2 Yoiar) + ps 






v4 = 


PM = f 4a;^(i-c)^^(1°sW - 1) + ^1?-^ + Ai2log(r) + ^3 






^ < 


Pm = Pinhom + Pi Ia{-ar) + p2 ■ A'o (ar) + ^3 


n = 3; 


flock 


^ > 


Pf = Pi sin(ar)i + p2 






^ = 


PF = pir^ + P2 






^ < 


Pf = Pi sinh(ar)i + p2 



with A = P 



\A\, and p satisfying p > 0, J pdx = 1. 



Proof. Let us define the operators £1 := A — fc^/, £2 A — | 
both operators to the equation and obtain 

.1 



C2CiiW*p) = {C2CiW)*p^ X{-CCiC2V 
= A ( -CP-^AS 



We apply 
C2CiVir)]*p 



Ck'F 



= A(l - Cr-^)Ap + A l^Ck'r-' -—jp = C2C1S 
using Lemma [TJ Hence, p should satisfy the following equation in its support: 



S + A6 - —5 ] * p 



k"^ 



Ap±a^p = - 



1 



1 



A 1 - 



-C2C1S, 



(14) 



with 



1^1 and 



A 



Ck^F- 



1 



1 - p _ (Jln ' 

resulting in the Hclmholtz equation for A > 0, the screened Poisson equation 
for A < and the Poisson equation for A = with radially symmetric inho- 
mogeneous right-hand side. Therefore, the solution to p4)) writes as a general 
solution of the homogeneous problem given by a linear combination of the fun- 
damental system in Remark [2] plus a particular solution of the inhomogeneous 
problem. 

The right-hand side of depends on the type of solution we wish to com- 
pute. For flocks in any dimension, s{r) is a constant function, and then we have 
\{i-cir'-'i) ^2-Cis{r) = D. Therefore, the inhomogeneous solution of (fT4)) for 
A ^ with unknown constant right-hand side D is 



.D 

PinhoTa^Ay) l^supp p ) 
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For mills and n = 2 we have s{r) = D + ^ log(r) to obtain 

■ log(r) + D (15) 



^ C2C1 



A(l-C) 



D + ^ log(r) 



/3 

fc4 



since log(r) is the fundamental solution of the Laplacian and its Dirac delta 
disappears and we look for mill solutions on an annulus (see Section [21 ©)■ 
Therefore, the inhomogeneous solution of ([T?)) with right-hand side can 
also be written explicitly. Again since log(r) is a fundamental solution of the 
Laplacian and the support of the solution is assumed not to contain the origin, 
it states 



a, , ^ D 



Pinhom,A(?') = , 212(^ ~~ log(^) + — SUpp p for A 7^ 0. 



Finally, in case A = 0, the inhomogeneous solution for the flock case is 

Pinhom,0 1 1 9 r, ' 

[g£)H ,n = 3 

whereas in the mill case it reads 

by using the fundamental solution of the Laplace operator. Putting together 
the inhomogeneous solution with the homogenous part leads to the claim of the 
theorem. For flocks, the space of candidate solutions is of lower dimension, since 
singularities at the origin are excluded. 

The coefficients {^1,^2) or {ni,^2,f^3) have to be computed numerically under 
the constraint that the solution has to be non-negative, has to contain unit 
mass, and has to solve the original equation (|13p . but only on its own support 
which is a priori unknown. To achieve this, we now need only to evaluate the 
convolution integral in (|13|) in a constrained optimization method rather than 
its inversion. 



Remark 3. Finally, we show that the radius of the support R and the constants 
(111,(12) o,fe connected by an explicit nonlinear identity in the particular case of 
3D flocks. By plugging the definition of the Quasi-Morse potential in 3D \12\ 
into ()lip . then 
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By looking up in the table of Theorem[Jl we have the explicit expression of ppir) 
for flocks with A> and n = 3. Straightforward computations lead to 

j\{r,s)pp{s)ds-^{Cl^ 



where A ~ was used, and with 

A(C, /) = fiik'^al cos Ra + p,2ka^l^R + l^fi2a^ + H2lk'^ + /iifc^ sin Ra + p,2k^R- 

Therefore, the existence of a flock solution is equivalent to the conditions A(C, I) = 
and A(l, 1) = 0, or equivalently 

fk'^al cos Ra + k^ sin Ra ka^PR + a^P + Ik"^ + k^R\ //iA _ /0\ 
\^ fc^a cos Ra + k^ sin Ra ka^R + + fc^ + k'^R j \/i2/ ~ ^0/ ' 

A necessary condition (pp can still he negative) is that the determinant of the 
matrix on the left hand side is zero, i.e., the existence of the solution of the 
nonlinear equation for R 

_a k^R- a^{P + l + klR) 
tani?a - -^^2^^^2(/2 + ; + l) + p- 

So that, for flock solutions in 3D we only need to check for radii verifying this 
last identity. 

In the next section, we will show an algorithm to solve this problem and present 
the numerical results. 



-1) = 



XCl^ 
A 

~ P(fc2 + a2) 



-A{C,l)e 



A(l,l)e 



_.j^ sinh(fr) 
r 

sinh kr 



r 



4. Numerical investigations 

Theorem [1] shows that solutions of ([T^ arc solutions of with the con- 
straints of positivity, unit mass, and compact support on an annulus. Therefore, 
we now propose an algorithm that numerically determines the support and lin- 
ear factors Pi of the stationary flock and mill solution. We will also present 
results which arc compared to particle simulations. 

4.1. The algorithm 

Let parameters n,C,l,k,a, P be fixed and phom denote the homogeneous 
solution dependent on dimension as in Theorem [T] In search for the support of 
the solution, we set the parameter Rmax as an upper boundary on the support 
size the algorithms shall consider. We can ensure that this is no restriction to 
the final result by setting Rmax large compared to the characteristic shape of the 
potential. Furthermore, denote Ar a discretization parameter and {ro, . . . , r^} 
an equidistant discretization of a chosen support supp(p), with ri+i — = 
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Ar. Numerical approximations of functions F(r) on the discrete radial grid are 
denoted with F. Our first algorithm determines the best possible solution for 
one particularly chosen support B(Ri,Rr). We aim to find linear coeflacients 
(/L<i,yU2) (or {fii, IJ.2, l^s) respectively), which solve the integral equation the 
best possible way. Non-negativity and unit mass of p are hard constraints, 
whereas the deviation W*p — s serves as the objective function the coefficients 
shall minimize. 

Algorithm 1 (for flocks). 
Input : fixed support B{0,Rr) 

- For convolving functions with the potential, compute a matrix H s.t. 



W -k p = Hp according to Section 12.21 

- Evaluate the convolution of the basis functions phom and 1 on supp p : 

:= Hphom,g'^ := Hi. 

- To fit the right hand side s(r) = on the support, we chose coefficients 
such that s(r) ~ D at the two end points ri, rpf. That is, solving 

gj, g%)^-''^[l 
setting D ^ I temporarily. 

- By linearity of H, we set p :— p-(/iconst.iPiiom + /^const.2) with M normalizing total mass. 

- Since we have only ensured (jl3p to hold at two points, we measure deviation of 
Hp from s(r) (here, an arbitrary constant) on the whole support as 



1 

Output : e, p, s if p > 0, error message if p ^ 



Hp-^ I Hpdf 



dr. 



For the case of mills, we proceed analogously, but we have to take into 
account the fixed inhomogeneous solution and three basis functions. 
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Algorithm 1 (for mills). 
Input : fixed support B{R„i, Rm) 

- For convolving functions with the potential, compute a matrix H s.t. 



W -k p = Hp according to Section [ 

- Evaluate the convolution of the fixed inhomogencous part /Oinhom.A on supp p 

and set Sinhom := -ffpinhom.A- 

- Define the remainder of the right-hand side as Sicm s — Sinhom, 
which has to be fitted by the convolution of the basis functions. 

- To do so, evaluate Ja{ar),YQ{ar) and 1 on suppp: 

■.= HJo,g^ :=HYo,g^ := HI. 

- Giving three basis functions, we pick three points ri,rj with j = \_N/2\,rN 
and interpolate both the remainder Sj-em and the free constant, which is 
temporarily set to 1. We solve 



Mconst 

- By linearity of H, we set prcm := ^h■cm.lJo + Mrcm,2>o + Mrcm,3 and 

Pconst := A*const,l>/o + Mconst,2^0 + /^const,3- 

- Our last degree of freedom is the free constant on the right hand side s(r), 
which we use to normalise mass. The candidate density is 

- - I - I - -.1 1 — ™(Prcm) — m(Pi„hom ,4 ) 

P /Oinhom.A + Prom + 7Pconst With 7 := m(pcon.t) ' 

- We penalize deviation of Hp from s(r) on the entire support as 

1 



9l 


9l 


9f 
9^ 


9] 


9] 


9n 


9n 


9n 








9\ 


9l 


9l 


9] 


9] 


9^ 


9n 


9% 


9n 




ei 



Hp- s - — — I {Hp - s)dr 



dr. 



Rm — Rn 

Second, since s(r) is concave, we penalize numerical convexity of s by 



62 := J X[s">o] sdr 

The total penalty value is the sum of ei, 62- 
Output : e = ei + 62, p, s if p > 



Now, we search the minimizer of the error function e over a test set of supports, 
given by the pre-defined discretization Ariand maximal support size. Repeating 
Algorithm [T] over the set of test supports provides a minimizer of the penalty 
function. For flocks, the number of tested supports is « , for mills w 

i ( ^Ar-T ) ■ snhance the speed of numerical computation, we first compute a 
solution based on a coarser discretization length Ar2 = mA^j for some integer 
m. Then, the obtained minimizer is used as the center of a local refinement 
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search with a fine diseretization length, as illustrated in Algorithm [21 
Algorithm 2. 

- Choose a coarse grid size Ar2 such that an iteration of Algorithm [T] over all 
test supports is reasonably fast and, as a solution, obtain the support -6(0, Rp) 

(or B{Rm,RM)) for mills. 

- Vary this support locally up to a fixed parameter c with a fine discretization 
Ari ^ Ar2, re-run Algorithm [T] restricted on \Rp — Rp\ < c 

{\Rm - Rni\ < c, \Rm - Rm\ < c for mills). 

- Obtain the minimising results supp(p), p, s and e. 



Naturally, the matrix H is not recomputed in every iteration but constructed 
once for the largest support and inherited. The choice to fix a functional equal- 
ity on the points which are most left, most right and for mills central on the 
chosen support is arbitrary. We say that no compact solutions are found in our 
computations, if our algorithms deliver i?rnax as the error minimizer, no matter 
of its value. The convergence of the algorithm for Ar — >■ if compact solutions 
are found will be demonstrated together with the results of the next subsection. 

4.2. Flocks in 2D 

We start our presentation of numerical results with the aligned fiock in two 
dimensions. Our standard example is the configuration C = Z = 0.75, k = ^ 
as in Fig. [TJ The stationary aligned flock state is independent of A, a, /?, yet 
emergence of flocks in particle simulations depends on these parameters and 
suitable initial conditions. An exemplary convenient choice is a = 1, /3 = 5, A G 
{100, 1000}. The observed flock of aligned particles is illustrated in Fig [2^ for 
N = 400 particles. In Fig. [^b, the result of our investigations is compared to 
the empirical radial density obtained from a particle simulation with N ~ 30000 
agents. The empirical radial density is obtained by collecting particles in radial 
bins and dividing by the Jacobian of the radial transformation. We see that the 
continuous solution matches the particle density and convergence is expected as 
iV — oo. While the numerical cost of full particle simulations is at least 0{N'^), 
the computational effort of the presented method scales quadratically with Ar, 
as illustrated in Fig. [3}d. In Fig. [3^ we show the convergence of our algorithms 
as Ar — >■ 0. One observes that the support is estimated well for coarse grid sizes, 
whereas the correct radial density is established with flner discretizations. The 
minimal error values of Algorithms [1] [2] are listed in Fig. [Sb- The advantages 
of the presented solution are continuity, dramatic reduction of the numerical 
cost, fast convergence, and an explicit expression of the radial density as, in this 
example, a combination of Bessel's J-function and a constant. 

Concerning the potential parameters, the area of relevant short-term re- 
pulsion and long-range attraction shapes divides into two subregions based 
on the results of section [31 as illustrated in Fig. [H In region I with C > 
1,1 < 1,CP < 1, the potential is catastrophic, ^ > (from Theorem [T|) and 
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(a) Flock emerged in a particle simulation with (b) Radial flock density: Continuous result vs. 
N = 400 particles empirical density (A'^ = 30000 particles) 

Figure 2: Two-dimensional aligned flocks emerge for the Quasi-Morse poten- 
tial. The resulting continuous radial density of Algorithms I1I2I matches the 
empirical distribution obtained from particle simulations. The stationary flock 
has the form pp = /ii Jo(ar) + ^2 with, in this case, fii sa 0.2356, /i2 ~ 
0.018,^ = 1.5, Rp ~ 1.31 (Quasi-Morsc potential parameters in use are 



(1) 


A 


=0.1 


(2) 


A 


=0.05 


(3) 


A 


=0.01 


(4) 


A 


=0.0025 



Ar 


error e 


computation time 


0.1 


3.54e-05 


0.76s 


0.05 


1.36e-05 


2.85s 


0.01 


3.99e-06 


69.1s 


0.0025 


9.97e-07 


1125s 



(a) Continuous solution pp for varying Ar (b) Minimal error value and computation times 

Figure 3: Algorithms II 121 converge as Ar — > if a compactly supported flock so- 
lution exists. Four resulting densities are shown for Ar £ {0.1, 0.05, 0.01, 0.0025} 
together with the minimal error value of the algorithm and the corresponding 
computation time. 
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Quasi-Morse potential has no minimum 




\ Region II: A <0, 

\ No continuous compactly supported 
liocks, no mill solutions 






Separatrix: 








as region II 






Region 1: A>0, 

Compactly supported flocks and mill solutions 



Figure 4: Numerical phase diagram of the Quasi-Morse potential in 2D: The 
biologically relevant scenarios decompose into two subregions. Region I: A > 
0, continuous compactly supported flocks. Region II: A < 0, no compactly 
supported continuous solutions, flocks only emerge on particle level. The same 
division of regions applies to the mill solutions. 

compactly supported continuous flock solutions are found. In region II with 
C > 1,1 < 1,CP < 1,A < (from Theorem [ij, and no compactly supported 
solutions can be found. No solution of the algorithm indicates H-stability since 
in this case particle simulations do show flocks whose support diverges when 
N ^ oo. The presented method faces numerical difficulties for catastrophic 
potentials A > with CP ~ 1, where it eventually breaks down not converging 
to a compactly supported flock. Similarly, particle simulations are not fully 
reliable in this limiting cases. However, thanks to our computation in Section 
Owe are able to consider the exact separatrix case CP = 1,C > 1,1 < 1,A = 0: 
Here, no compact solutions arc found. Our numerical findings are illustrated in 
Fig. |4l We emphasize that based on the reported simulations, wc conjecture 
that compactly supported flock solutions exist only in the catastrophic regime, 
^ > 0. 

4.3. Mills m 2D 

The Quasi-Morse potential is able to produce rotating mill states in particle 
simulations, just as the original Morse potential. We choose the same config- 
uration as in Section 14.21 with A = 100 and show the mill emerging from a 
particle simulation in Figure [5ji. The resulting mill solution of our algorithms 
is illustrated in Figure [2}d, together with a comparison to an empirical density 
from a particle mill with N = 16000 agents. Again, our result is confirmed by 
the particle simulation and support as well as the density shape agree perfectly. 
The stationary rotating mill is a weighted sum of Bessel's J and Y functions, 
the inhomogeneity phom and a constant. The convergence of Algorithms [U [2] in 
the mill case is shown in Figure El As for flocks, the computational costs are 
minimal compared to a full particle simulation. For the existence of compactly 
supported mill solutions, the parameter diagram on Figure |4] applies just as for 
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(a) Mill emerged in a particle simulation with (b) Radial mill density: Continuous result vs. 
N = 400 particles empirical measure (A'' = 16000 particles) 



Figure 5: Rotating mills emerge for the Quasi-Morse potential. As for flocks, 
the resulting radial density of Algorithms I1I2I matches the empirical distri- 
bution obtained from particle simulations. The mill solution has the form 
Pm = Pinhom,A + Mi Jo{ar) + fi2 ^o(a»') + M3 with, in this case, w 0.1708, fi2 ~ 
0.0468, = 0.0320, A = 1.5,supp^^^ w 5(0.47,1.57) (Quasi-Morse potential 
parameters in use are C = ^,1 = 0.75, k = ^, others are a = 1, /3 = 5, A = 100). 




Figure 6: Algorithms I1I2I converge for the mill case as Ar — > 0. 
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(a) mill solutions for varying ^ (b) support of mill solutions for varying ^ 

Figure 7: Quasi-Morse potentials with identical shape parameters C,l,k result 
in mill solutions with different support sizes and densities, depending on the 
ratio of potential factor and squared stationary speed of the mill. 

flocks. In region I, continuous solutions can be found, whereas in region II and 
the separatrix CP = 1 no such mills can be found. In particle simulations, 
we there see either a crystal-like arrangements or "finite particle" flocks as in 
Section |4?2] Next we study the impact of parameters a,/3, A on the stationary 
mill solution, which enter the solution solely in the joint quotient Hence, 
for a potential multiplied by a factor A, the mill solution will stay the same, if 
the preferred speed of particles is multiplied by by any suitable change of a 
and/or /?. In Figure[7^, we show several mill densities for our standard potential 
configuration and ^ S {350, 500, 1250, 2500, 5000, 12500}. The support of mill 
solutions is plotted against ^ in Figure [TJj. 

4 ■4- Flocks in 3D 

The introduction of Quasi-Morse potentials enables us also to study flocks in 
three space dimensions. As we have mentioned in Section |3l the area of admissi- 
ble parameter configurations is smaller than in the 2D case, as illustrated in the 
parameter diagram. For our example, we set C = 1.255, Z = 0.8, fc = 0.2, A = 
5.585 and plot the resulting potential shape in Figure |H^. A three-dimensional 
flock resulting from a particle simulation is shown in Figure |8}d. With the help 
of Algorithm[T]the continuous radial flock density is computed as a linear combi- 
nation of SiSiii and a constant. Notice that due to Remark|31 Algorithm [5] is not 
needed. Also in three dimensions, the empirical density of a particle simulation 
matches our result, as illustrated in Figure[51:. Concerning the existence of flock 
solutions in dependence of the shape parameters C and /, we get an equivalent 
picture as in two dimensions (see Figure[8jl): Though different in shape, the area 
of biologically relevant shapes is divided into two subregions by the separatrix 
Cl^ 1. In region I, continuous compactly supported three-dimensional flocks 
are found, not so in region II, which again indicates H-stability. Here, flocks do 
appear but their support increases with the total number of agents N . In the 
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special case of the separatrix, which can be investigated with the computation 
of the case A = in Section [3l no flock solutions are found. 

5. Discussion 

Quasi-Morsc potentials fulfill three properties desirable from biological mod- 
eling: short-term repulsion, long-term attraction and vanishing interaction at 
infinity. Using Quasi-Morse potentials instead of the standard Morse potential 
makes, in our view, hardly any difference in terms of biological modeling. The 
stronger singularity at the origin for n ^ 2 might even be desirable in order to 
enforce repulsion. Though the special functions involved for n = 2 may seem not 
as convenient to work with as the exponential function, the existence of contin- 
uous, compactly supported stationary states itself make Quasi-Morse potentials 
a good choice for further studies of the models discussed in the above. Our 
results are, to the best of our knowledge, one of the first of its kind for explicit 
solutions of flock and mill patterns in two or three dimensions. The strategy 
of building up potentials from solutions of certain partial differential equations 
might work in other cases as well and form one tool in the effort to understand 
the equilibria of interaction potentials. However, the techniques applied here 
are of no help for general potentials, such as classical Morse. With a variety of 
potentials suggested (see the discussion in Section [2]), the problem of choosing 
the best suited one for a particular biological application becomes increasingly 
evident and should be a topic of future research. 

6. Conclusions 

In this paper we have introduced the Quasi-Morse interaction potentials for 
a second-order model of self-propelled interactive particles. The Quasi-Morse 
potentials lead to the emergence of flocks and mills, similar to the standard 
Morse potential. We have shown that the radial densities of these stationary 
states are (afiine) linear combinations of two or three elementary functions, 
which are chosen with the respect to the three subcases A > (catastrophic), 
A = (separatrix) oi A < 0. In order to determine the correct scalar coefficients 
and the a priori unknown support, we have developed a numerical algorithm that 
does not use time evolutions in Section 21 We have illustrated our result with 
examples for flocks and mills in two dimension and flocks in 3D. In all cases, our 
findings are convincingly verified by corresponding particle simulations. With 
our algorithm, we find that for all coherent patterns, only the catastrophic 
scenarios A > lead to continuous compactly supported solutions. 

7. Appendix 

This appendix is devoted to the proof of Corollary [2j Denote by Vk (r) the 
potential in ([T2|) for a given fc > 0. Notice that Vfc(r) = fc"~^Vi(fcr), hence we 
set fc = 1 without loss of generality and, from now on, wc drop the index k for 
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(a) Quasi-Morse potential in 3D 



(b) 3D flock emerged in a particle simulation 
with N = 200 particles 




(c) 3D radial flock density: continuous result vs. (d) Parameter diagram 

empirical measure (TV = 35000 particles) 

Figure 8: The Quasi-Morse potential in three dimension is able to produce 
aligned flock sohitions. The continuous radial density can be expressed as pF = 
/xi-sin(ar)i+/i2-l with, in this case, fii ~ 0.3574, /i2 « 0.0052, « 0.725,^ = 
5.585. Our result is verified by comparing to the empirical density obtained from 
a particle simulation, (a) Exemplary potential shape, (b) Flock emerged from 
3D particle simulation, (c) Continuous flock solution vs empirical measure, (d) 
Parameter diagram of biological relevant configurations 
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simplicity. 

Let us first show the assertion about the unique minimum of the potential. 
This property is desirable from the biological point of view to set a typical length 
scale for the distance between agents. 

Let U{r) = T/(r)-Cy(f), then U'{r) = V (r) - {^) , a.nd the necessary 
condition for a local extremum can be stated as finding r such that 

Straightforward computations lead to 

- nrwm = j^^^{iv'{r)v{^j) - rim-)), 

where we used that y is a solution of ^^^^V^'(r) + V"{r) — V{r). 

We next check that log(— F(r)) is a convex function of r for n = 1,2,3. 
Indeed, if n = 1, then log(— F(r)) is affine; if ri = 3, then log(— V^(r)) = — logr — 
r — log(47r). In the case n = 2, we have 

(log(-F(r)))" = K2{x)lMx) - 2K,{xf 



2A'o(.t)2 

and we can use the inequality Kq{x)^ + K2{x)Kq{x) — 2A'i(x-)^ > (which can 
be verified numerically). 

Thus, if / < 1, we have > Since V{r) < 0, V'{r) > 0, this 

implies V'{^)V{r) -V'{r)Vi^)>0, and therefore 



Similarly, if Z > 1, we obtain h'{T) > 0. 
Further, it is directly checked that 

lim h(r) = Cr-^ and lim h(r) = I ° !^ ! ^ ! ■ 

r^0+ ^ ' r^oo ^ ' \ +0O if / > 1 

Thus, the equation hir) = 1 has no solution in the cases Cl"^^ < 1, Z < 1 or 
(jin-2 > 1 I > 1 and a unique positive solution in the cases Cl^^^ < 1, Z > 1 or 
CZ"-2 > 1, ? < 1. Recalling that U'{r) = V'{r){l - /i(r)), we see that of the last 
two cases, the former corresponds to a local maximum of U{r) and the latter to 
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a local minimum. Since the minimum is unique, then it is a global minimum. 
Concerning the second assertion, by construction we have 




y(|x|)d.i; = — 1 for all n , 



since V is the fundamental solution of Au — u — 6q as stated in Definition [31 
Therefore, we get 

/ U{\x\)dx= [ V{\x\)-CV{\x\/l)dx = -l + Cl'', 

JR" JR" 

which is negative for CI" < 1, and thus U is catastrophic (see [s^, p. 37). 
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